Introduction

Bell pepper (Capsicum annuum L.) is an annual and herbaceous plant that belongs to the family Solanaceae such as tomato and potatoes, and is one of the most economically important crops in the world. In 2017, bell pepper was considered the third vegetable with the highest production worldwide, with an estimated contribution of 36 million tons. Since Capsicum grows in tropical and even temperate regions, diverse abiotic stresses, such as salinity, temperature, drought, flood, UV radiation and heavy metals, may affect its growth, causing 50 to 70% yield losses worldwide (Chugh et al. 2018).

In bell pepper, the temperature greatly affects its production, in which the optimal temperature ranges from 21 to 27ºC, while lower temperatures affect its growth and reproduction (Pressman et al. 2006). Several studies have shown that cold induces numerous morphological, biochemical and molecular changes in C. annuum. Mercado et al. (1997) observed a decrease in height, number of leaves and leaf area, while the content of carbohydrates and soluble proteins were increased. In leaves, exposure to 8ºC increases the levels of antioxidant compounds as ascorbate, glutathione and NADPH-generating dehydrogenases (Airaki et al. 2012). Likewise, Guo et al. (2012) showed that cold (10/6°C) increased H₂O₂ and malondialdehyde, indicating cell membrane damage, which consequently triggers an increase of enzymatic activity of glutathione reductase, dehydroascorbate reductase, monoDHAR, guaiacol peroxidase and ascorbate peroxidase. In pepper seedlings, cold treatment increased the accumulation of total soluble proteins, proline and phenolic compounds in stems, while decreased the content of chlorophyll (Koç et al. 2010). Molecularly, several transcription factors are induced upon exposure to cold stress, including EREBP (CaEREBP-C1 to C4), WRKY (CaWRKY1), bZIP (CaBZ1) (Hwang et al. 2005), NAM, ATAF1/ 2, CUC2 (Hou et al. 2020) and ERF/AP2-type (CaPF1) (Yi et al. 2004), in which heterologous overexpression of CaPF1 increased tolerance against freezing and resistance to pathogens in Arabidopsis (Yi et al. 2004), while overexpression of CaNAC064 increased tolerance to cold stress (Hou et al. 2020).

On the other hand, ultraviolet-B radiation (UV-B), corresponding to the high energy (280–320 nm) of daylight, has a great impact on plants. In bell pepper leaves, UV-B was found to increase proline, quercetin, rutin and anthocyanin, while the content of chlorophylls and carotenoids were reduced (Mahdavian et al. 2008). Moreover, Rodríguez-Calzada et al. (2019) reported an increased expression of the phenylalanine ammonia lyase (PAL) and chalcone synthase (CHS) genes, related to the accumulation of chlorogenic acid, luteolin 8-C-hexoside in response to UV-B. Another study, Lai et al. (2011) identified 183 differential expression genes, related to carbohydrate metabolic process, protein modification process, catabolic process and cellular homeostasis.

In nature, the combination of two or more stresses is common, and plant responses induced by combined stressors are largely controlled by cross-talk between different sensors and signal transduction pathways, which can activate or inhibit each other (Mittler and Blumwald 2010; Atkinson and Urwin 2012; Suzuki et al. 2014). Despite the advances in understanding the molecular regulation in UV-B or cold stress, a few studies have been conducted to assess the combined effect of these abiotic factors in plant stress responses. In this regard, León-Chan et al. (2017) showed that UV-B and cold induced the degradation of chlorophyll and accumulation of carotenoids, chlorogenic acid, apigenin and luteolin glucosides in comparison to each abiotic stress. Further, transcriptional analysis showed the upregulation of flavanone 3-hydroxylase (F3H) gene indicating the activation of flavonoid biosynthetic pathway in response to UV-B and cold in bell pepper stems, while flavonoid-3', 5'-hydroxylase (F3´5´H), dihydroflavonol-4-reductase (DFR) and anthocyanidin synthase (ANS) were more strongly induced separately in UV-B or cold treatments (León-Chan et al. 2020). Nonetheless, changes in global gene expression patterns in response to combined UV-B and cold is relatively unknown. In an attempt to gain deeper insights about the temporal dynamic of genes and pathways modulated by these combined stressors, we analyzed transcriptional changes using the RNA-seq analysis to provide relevant information about the genes and the pathways that participate in stress-induced responses. Hence, the aim of this study was to analyze the transcriptomic profile of C. annuum stems in response to combined UV-B radiation and cold stress at different times, to provide new insights about the specific genes and pathways involved at early, intermediate and late plant responses.

Materials and Methods

Plant material and growth conditions: Commercial bell pepper seeds Canon cv. (Zeraim Gedera Syngenta; Israel) were germinated and maintained as previously described (León-Chan et al. 2017). Twenty-eight days after sowing (DAS), bell pepper plants were put into a plant growth chamber (GC-300TLH, JEIO TECH; South Korea) at control conditions, which consisted of a 12 h photoperiod (from 6:00 to 18:00 h) of PAR radiation (972 μmolm^-2s^-1), temperature of 25/20°C (day/night) and relative humidity of 65% for three days. For treatment of UV-B and cold, temperature was adjusted at 15/10°C the previous night (day 30 at 18:00 h) and Capsicum plants were irradiated with PAR for 6 h (from 06:00 to 10:00 and 16:00 to 18:00 h) and UV-B irradiation (72 kJ·m²) for 6 h (from 10:00 to 16:00 h), and this was maintained until sampling (day 31 and 32). For sampling, stems from 10 bell pepper plants were collected at 0, 1, 3 and 25 h after stress exposure by duplicate, frozen in liquid nitrogen and stored at -80°C.

Total RNA isolation and library preparation: Treated and control plant stems were collected and subjected to total RNA isolation. Stems were pulverized with liquid nitrogen and total RNA was isolated from 50–100 mg of tissue with Trizol reagent (Ambion, Life Technologies, U.S.A.) according to the manufacturer ́s instructions with the following modifications: for precipitation step, we replaced 0.5 mL of isopropyl alcohol, with a mixture of 0.25 mL of isopropyl alcohol and 0.25 mL of 7.5 M lithium chloride; finally, RNA washes with 75% ethyl alcohol were carried out twice. Genomic DNA was removed with Turbo DNA free kit (Invitrogen, Life Technologies, U.S.A.). RNA concentration was determined using NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, U.S.A.) and RNA integrity was analyzed by agarose gel electrophoresis. RNA of acceptable purity and integrity (A260/A280: ≥ 1.8; RIN ≥ 8) was used to prepare cDNA libraries of 150 paired-end readings in the Illumina TruSeq library system. The concentration of two libraries was determined by fluorometry at Qubit (Life Technologies). Later, the libraries were sequenced on the Illumina NextSeq-500 platform according to the sequencing service provider, National Laboratory of Genomics for Biodiversity (LANGEBIO) Unit CINVESTAV-IPN; Irapuato, Guanajuato, Mexico.

Data processing and DEG identification: The quality of raw reads was visualized using FASTQC program, and then trimmed using Trimmomatic with the following parameters: quality score of 30 (SLIDINGWINDOW:4:30) and minimum reading length of 20 (MINLEN: 20). Afterwards, the trimmed reads were aligned to the pepper reference genome (Pepper Zunla 1 Ref_v1.0, https://www.ncbi.nlm.nih.gov/genome/?term=txid4072[orgn]) using HiSAT2. Gene expression levels were calculated by counting the number of mapped reads per annotated gene model using HTSeq-count, and raw read counts were normalized for RPKM (Love et al. 2014). For downstream analyses, differentially expressed genes (DEG) were determined using DESeq2 in R software (Anders and Huber 2010), where DEGs were considered with ≧ 1.5-fold expression with respect to the control and adjusted P value α ≦ 0.05. The Volcano plots, Venn diagrams and Cluster analysis were realized using pheatmap, EnhancedVolcano and VennDiagram package in R software (version 1.2.5001; http://www.r-project.org/).

Gene ontology and KEGG enrichment analysis: The GO enrichment of DEGs was performed in UNIPROTKB (https://www.uniprot.org/uploadlists/) and AgriGO (http://bioinfo.cau.edu.cn/agriGO/analysis.php) web-based tool for GO analysis. GO terms were performed with FDR ≦ 0.05. We carried out the statistical enrichment of the differential expression genes in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (α ≦ 0.05).

Results

Data processing and DEG identification: A total of 97,291,544 paired-end raw reads were obtained in this study. The quality assessment using FastQC showed an average 24,663,149 reads with a length of 150 pb and an average content of 51% GC, per sample. All raw reads from samples had quality levels with a Phred value between 14 and 36. After filtering with Trimmomatic, the samples were left with filtered reads with length > 20 and a Phred value ≥ 30 (Q30), preserving on average 39% of the total raw reads (Table 1). Read alignments had a mapping rate of 35.34 to 50.87% of total filtered reads.

For the combined treatment of UV-B and cold at 1 h, 281 differentially expressed genes (DEG) were identified, of which 154 were up-regulated and 127 down-regulated (Fig. 1A); for 3 h, 280 DEG, of which 167 were up-regulated and 113 down-regulated (Fig. 1B); and for 25 h, 326 DEG, of which 138 were up-regulated and 188 down-regulated (Fig. 1C).

The Venn diagrams revealed that the gene expression profile differed significantly along the three treatments, showing that 29, 54 and 32 genes were up-regulated exclusively at 1, 3 and 25 h, respectively, and 66 genes were induced at all-time points (Fig. 2A). For down-regulated genes 29, 28 and 90 were exclusively observed at 1, 3 and 25 h after treatment exposure, in addition to 65 genes observed at all times of sampling (Fig. 2B). Interestingly, the 66 up-regulated genes present at all times, included genes such as APRR1, APRR5, chalcone synthase-1B, chalcone synthase-2 and chalcone synthase-J related to photoperiod and flavonoid synthesis; whereas the 65 down-regulated genes expressed at the different times of sampling were involved in diterpenoid, sesquiterpenoid and triterpenoid biosynthesis and defense against pathogens, such as beta-amyrin synthase, (-)-

Table 1: Statistics of raw reads filtering

Sequences	Total raw reads	Total filtered reads (Trimmomatic)	Q30 (%)	GC (%)
Ctrl A	22387045	8440660	38	51%
Ctrl B	23517336	9468927	40	51%
Treat 1A	24488878	9073092	37	51%
Treat 1B	26898285	10541133	39	49%
Treat 3A	24884652	9742863	39	51%
Treat 3B	20673273	7516910	36	49%
Treat 25A	30987008	13585897	44	42%
Treat 25B	23468715	9219918	39	51%

Ctrl A, Ctrl B: stem of control samples; Treat 1A, Treat 1B: stem exposed to 1 h; Treat 3A, Treat 3B: stem exposed to 3 h; Treat 25A, Treat 25B: stem exposed to 25 h. Ctrl: control; Treat: treatment; G-C: Guanine-cytosine

germacrene-D-synthase, cytochrome P450-82C4, PYL12, flower-specific defensin and zingipain. Cluster analysis revealed very different transcriptomic profiles underlying a marked differential gene expression at each time of sampling, including genes involved in plasma membranes, compound transport, chloroplast, cell wall, signaling and transduction of cellular signals, ROS oxidation, hormones and activity against pathogens (Fig. 3).

GO classification analysis of DEGs: For the combined treatment of UV-B and cold after 1 h of exposure, GO enrichment analysis showed that three categories for the cellular component, two for molecular function and three for the biological process, of which hydrolase activity (GO:0016787) and response to stress (GO:0006950) were statistically significant with 25 and 15 genes, respectively. For 3 h treatment, three categories were identified for the cellular component, two for molecular function and 14 for biological processes; from these, four categories were statistically significant: response to abiotic stimulus with 19 genes (GO:0009628), response to stress with 14 genes (GO:0006950), carbohydrate metabolic process with 8 genes (GO:0005975) and biosynthetic process with 18 genes (GO:0009058). For treatment at 25 h, two categories were identified for the cellular component, two for molecular function and five for biological process; from these, three categories were statistically significant: hydrolase activity with 36 genes (GO:0016787), response to stress with 17 genes (GO:0006950) and carbohydrate metabolic process with 18 genes (GO:0005975) (Fig. 4).

The response to stress (GO:0006950) category was significantly identified in all times, where genes were related to hormone biosynthesis, ROS oxidation and defense against pathogens, some genes are cytochrome P450 (98A2, CYP72A219 and CYP736A12), catalase, peroxidase, pathogenesis-related STH-2, RPP13 disease resistant and flower-specific defensin (Table 2). Genes grouped in hydrolase activity (GO:0016787) were found at 1 h and

Table 2: Genes identified in the response to stress category

Gene ID	Name	1 h		3 h		25 h
Gene ID	Name	LFC	FDR	LFC	FDR	LFC	FDR
107854492	catalase	3.85	0.02	4.69	0.00	3.88	0.02
107856092	peroxidase 45-like	2.05	0.03	1.34	0.36	2.18	0.02
107871732	cryptochrome DASH, chloroplastic/mitochondrial	2.09	0.00	3.06	0.00	1.70	0.00
107844023	cytochrome P450 98A2-like	1.10	0.01	1.51	0.00	1.38	0.00
107878596	cytochrome P450 CYP72A219-like	1.54	0.00	2.11	0.00	1.03	0.08
107850965	cytochrome P450 CYP736A12-like	2.09	0.00	1.64	0.00	2.20	0.00
107863949	linolenate hydroperoxide lyase, chloroplastic	1.33	0.01	1.52	0.00	0.57	0.51
107877227	cytochrome P450 72A15-like	-4.84	0.01	-5.76	0.02	-0.39	0.93
107863881	cytochrome P450 82C4-like	-2.11	0.04	-3.19	0.00	-2.11	0.03
107870440	disease resistance protein RPP13-like	-1.48	0.07	-1.26	0.20	-1.65	0.03
107864567	pathogenesis-related protein STH-2-like	-1.50	0.00	-1.08	0.00	-0.95	0.01
107850294	kirola-like	-0.69	0.31	-1.01	0.07	-2.66	0.00
107877005	flower-specific defensin-like	-1.67	0.04	-2.08	0.01	-2.77	0.00
107863162	RNA polymerase sigma factor sigE, chloroplastic/mitochondrial	1.50	0.00	1.90	0.00	1.11	0.02
107875362	E3 ubiquitin-protein ligase CHIP	1.10	0.01	1.52	0.00	0.91	0.05
107865651	ethylene-responsive proteinase inhibitor 1-like	-1.25	0.55	-2.06	0.23	-2.52	0.05
107850595	dnaJ protein homolog	-0.91	0.03	-0.74	0.15	-1.72	0.00
107843192	protein ROS1-like	-1.09	0.03	-1.52	0.00	-0.99	0.06
107879996	Fanconi anemia group I protein	0.76	0.56	0.88	0.52	1.55	0.04
107864208	phosphate transporter PHO1	1.72	0.00	1.51	0.01	1.53	0.01
107848500	bidirectional sugar transporter N3-like	1.31	0.01	1.57	0.00	1.13	0.04
107845990	pyruvate decarboxylase 1	1.53	0.00	1.32	0.00	1.26	0.01
107859400	allantoinase	1.75	0.00	1.33	0.04	0.87	0.28

LFC: log2 fold changes

Fig. 1: Volcano graph, genes differentially expressed at 1 (A), 3 (B) and 25 h (C) in response to combined stress UV-B radiation and cold

25 h, corresponding to carboxylesterase 8, vicianin, zingipain, endochitinase, ABC transporter (B, C and G), acyl-thioesterase 1/2 and phospholipase D, which participates in defense against pathogens, plasma membranes and transport of compounds (Table 3). Besides, the carbohydrate metabolic process (GO:0005975) related to changes in the cell wall was important at 3 h and 25 h, finding genes such as β-D-xylosidase 2, β-galactosidase, pectinesterase, inositol oxygenase and endoglucanase (Table 4). Meanwhile, response to abiotic stimulus (GO:0009628) and biosynthetic process (GO:0009058) were only found at 3 h; interestingly, the genes identified in these two categories are related to photoreceptor activity, protection of chloroplasts and flavonoid biosynthesis, some genes were ultraviolet-B receptor UVR8, adagio 3, stress enhanced, dehydrin, sigma factor, chalcone synthase J, chalcone synthase-1B and chalcone synthase-2 (Table 5 and 6).

Fig. 2: Venn diagrams showing the shared differential number up-regulated (A) and down-regulated genes (B) between 1, 3 and 25 h

Fig. 3: Cluster analysis of differential genes at 1, 3 and 25 h after combined cold and UV-B treatment

KEGG analysis of DEGs: Regarding the relevant role of UV-B and cold in the modulation of metabolism revealed by GO enrichment, we decided to analyze DEG using KEGG enrichment map. Our analysis showed that DEG belonging to the circadian rhythm-plant and flavonoids biosynthesis were the most enriched among the 10 pathways identified to up-regulated genes at 1 h (Fig. 5A), while no pathway was significant for the down-regulated genes (Fig. 5B). In addition, the enriched pathways at 3 h of exposure to combined treatment primarily were circadian rhythm-plant and flavonoids biosynthesis for the up-regulated genes, both statistically significant (Fig. 6A), meanwhile for the down-

Table 3: Genes identified in the hydrolase activity category

Gene ID	Name	1 hour		3 hours		25 hours
Gene ID	Name	LFC	FDR	LFC	FDR	LFC	FDR
107862340	ABC transporter B family member 25-like	1.55	0.00	1.38	0.01	1.56	0.00
107855690	ABC transporter C family member 14	1.52	0.00	1.21	0.00	1.69	0.00
107878709	dynein light chain 1, cytoplasmic	2.42	0.01	2.28	0.01	1.99	0.03
107845681	myosin-2-like	-1.54	0.01	-1.56	0.01	-1.03	0.13
107839653	phospholipase D zeta 1-like	2.25	0.01	1.05	0.58	2.40	0.01
107843045	lipid phosphate phosphatase gamma, chloroplastic	1.53	0.01	1.37	0.03	1.38	0.02
107839030	acyl-protein thioesterase 1 homolog 1-like	1.63	0.03	1.79	0.02	2.03	0.00
107876040	acyl-protein thioesterase 2	1.75	0.00	1.83	0.00	1.46	0.00
107859400	allantoinase	1.75	0.00	1.33	0.04	0.87	0.28
107859674	ATP-dependent zinc metalloprotease FTSH 6, chloroplastic	3.81	0.00	3.77	0.00	3.21	0.00
107851284	fumarylacetoacetase-like	0.97	0.05	0.92	0.10	1.53	0.00
107869091	beta-amylase 3, chloroplastic	2.60	0.00	2.67	0.00	2.73	0.00
107878266	beta-amylase	1.57	0.05	1.25	0.22	1.46	0.07
107862757	alpha-galactosidase-like	-0.96	0.00	-0.79	0.04	-1.54	0.00
107867137	patatin-like protein 2	-1.84	0.00	-0.83	0.39	-1.20	0.07
107860321	PLP1; patatin-like protein 3	-1.80	0.13	-1.79	0.18	-2.26	0.02
107859333	phospholipase A1-IIgamma-like	N/A	N/A	-2.08	0.19	-2.48	0.03
107875477	ABC transporter G family member 31	-1.62	0.04	-0.74	0.65	-1.38	0.10
107870765	ABC transporter B family member 2-like	-0.37	0.85	-0.37	0.89	-2.32	0.00
107871378	pleiotropic drug resistance protein 2-like	-1.50	0.00	-1.40	0.00	-1.13	0.01
107872419	probable carboxylesterase 8	-1.95	0.00	-1.31	0.01	-2.50	0.00
107875683	CAF1; probable CCR4-associated factor 1 homolog 9	-0.96	0.30	-1.72	0.02	-1.82	0.00
107841124	basic 7S globulin-like	N/A	N/A	N/A	N/A	-3.36	0.04
107859803	acidic 27 kDa endochitinase	-1.33	0.02	-1.14	0.07	-1.50	0.00
107859806	basic endochitinase-like	-2.01	0.06	-1.72	0.17	-2.67	0.00
107856465	vicianin hydrolase-like	-1.00	0.07	-1.15	0.04	-1.82	0.00
107860257	zingipain-2-like	-1.67	0.02	-1.96	0.01	-2.59	0.00
107861184	zingipain-2-like	N/A	N/A	N/A	N/A	-3.20	0.03
107870929	serine carboxypeptidase-like 19	-1.19	0.22	-0.97	0.45	-1.63	0.03
107867007	subtilisin-like protease SBT1.2	-1.81	0.02	-0.64	0.73	-2.17	0.00
107840985	probable beta-D-xylosidase 2	-2.20	0.00	-1.86	0.00	-2.23	0.00
107854898	glucan endo-1,3-beta-glucosidase A-like	-3.16	0.00	-2.75	0.02	-2.14	0.05
107879143	glucan endo-1,3-beta-glucosidase, basic	-1.62	0.27	-1.51	0.38	-2.29	0.05
107840962	BG1; beta-galactosidase-like	-3.18	0.00	-2.09	0.00	-3.40	0.00
107861740	beta-galactosidase	-2.12	0.00	-1.56	0.00	-1.80	0.00
107863277	pectin acetylesterase 9	-1.78	0.01	-1.33	0.15	-0.90	0.37
107859553	pectinesterase-like	-2.07	0.00	-1.80	0.00	-2.20	0.00
107864477	ccel1; endoglucanase 18-like	2.85	0.04	2.46	0.13	3.02	0.02
107843046	DEAD-box ATP-dependent RNA helicase 57	2.07	0.00	1.62	0.03	1.77	0.01
107862137	nudix hydrolase 18, mitochondrial-like	1.51	0.00	1.02	0.13	1.37	0.01
107859272	xylem cysteine proteinase 2-like	1.66	0.03	1.60	0.05	0.88	0.43
107874054	probable ribonuclease P/MRP protein subunit POP5	1.61	0.02	1.24	0.18	1.53	0.03

N/A: these genes are not differentially expressed; LFC: log2 fold changes

Table 4: Genes identified in the metabolic carbohydrate process category

Gene ID	Name	1 h		3 h		25 h
Gene ID	Name	LFC	FDR	LFC	FDR	LFC	FDR
107869091	beta-amylase 3, chloroplastic	2.60	0.00	2.67	0.00	2.73	0.00
107862757	alpha-galactosidase-like	-0.96	0.00	-0.79	0.04	-1.54	0.00
107859803	acidic 27 kDa endochitinase	-1.33	0.02	-1.14	0.07	-1.50	0.00
107859806	basic endochitinase-like	-2.01	0.06	-1.72	0.17	-2.67	0.00
107859802	CAChi2; acidic endochitinase pcht28	-2.03	0.00	-1.54	0.01	-1.09	0.10
107856465	vicianin hydrolase-like	-1.00	0.07	-1.15	0.04	-1.82	0.00
107840985	probable beta-D-xylosidase 2	-2.20	0.00	-1.86	0.00	-2.23	0.00
107854898	glucan endo-1,3-beta-glucosidase A-like	-3.16	0.00	-2.75	0.02	-2.14	0.05
107879143	glucan endo-1,3-beta-glucosidase, basic	-1.62	0.27	-1.51	0.38	-2.29	0.05
107840962	BG1; beta-galactosidase-like	-3.18	0.00	-2.09	0.00	-3.40	0.00
107861740	beta-galactosidase	-2.12	0.00	-1.56	0.00	-1.80	0.00
107864477	ccel1; endoglucanase 18-like	2.85	0.04	2.46	0.13	3.02	0.02
107859553	pectinesterase-like	-2.07	0.00	-1.80	0.00	-2.20	0.00
107859925	GS; galactinol synthase 2	2.60	0.00	2.22	0.02	2.02	0.03
107850683	inositol-3-phosphate synthase	2.64	0.00	1.95	0.00	1.96	0.00
107840943	inositol oxygenase 4	-2.84	0.00	-3.70	0.00	-1.74	0.02
107867324	phosphoenolpyruvate carboxykinase [ATP]-like	-1.40	0.00	-1.45	0.00	-1.53	0.00
107878490	xyloglucan endotransglucosylase/hydrolase protein 15-like	-0.60	0.72	-1.07	0.38	-1.72	0.02
107860149	probable xyloglucan endotransglucosylase/hydrolase protein 7	-0.30	0.82	-0.52	0.60	-1.63	0.00
107847799	xyloglucan endotransglucosylase/hydrolase protein 31-like	1.22	0.01	1.07	0.06	2.53	0.00

LFC: log2 fold changes

Table 5: Genes identified in the response to abiotic stimulus category

Gene ID	Name	1 h		3 h		25 h
Gene ID	Name	LFC	FDR	LFC	FDR	LFC	FDR
107862948	(6-4) DNA photolyase	1.56	0.00	1.66	0.00	1.42	0.00
107838759	adagio protein 3	1.43	0.00	1.76	0.00	1.49	0.00
107851542	stress enhanced protein 2, chloroplastic	2.74	0.00	2.89	0.00	2.57	0.00
107871194	ultraviolet-B receptor UVR8	3.67	0.00	2.59	0.00	3.05	0.00
107873562	UV-B-induced protein At3g17800, chloroplastic-like	1.84	0.03	1.52	0.15	1.79	0.03
107842826	ultraviolet-B receptor UVR8-like	1.61	0.00	-0.36	0.86	1.46	0.01
107863294	low-temperature-induced 65 kDa protein-like	5.96	0.00	5.16	0.03	3.12	0.31
107860006	dehydrin HIRD12-like	1.89	0.06	2.08	0.03	0.71	0.71
107871210	dehydrin HIRD11-like	1.62	0.00	1.48	0.00	1.08	0.02
107858537	dehydrin Xero 1-like	2.27	0.00	1.82	0.00	0.91	0.25
107866811	Dhn; phosphoprotein ECPP44-like	1.54	0.00	1.50	0.00	0.80	0.01
107853534	mitogen-activated protein kinase kinase kinase ANP1-like	5.39	0.01	4.56	0.07	5.97	0.00
107855817	B-box zinc finger protein 32	3.41	0.28	5.52	0.02	2.25	0.59
107854515	protein PHYTOCHROME KINASE SUBSTRATE 4	-0.55	0.29	-1.62	0.00	-0.27	0.73
107862854	MKK1; mitogen-activated protein kinase kinase 9	-0.27	0.88	-1.88	0.01	-0.52	0.66
107863162	RNA polymerase sigma factor sigE, chloroplastic/mitochondrial	1.50	0.00	1.90	0.00	1.11	0.02
107848500	bidirectional sugar transporter N3-like	1.31	0.01	1.57	0.00	1.13	0.04
107875362	E3 ubiquitin-protein ligase CHIP	1.10	0.01	1.52	0.00	0.91	0.05

LFC: log2 fold changes

Table 6: Genes identified in the biosynthetic process category

Gene ID	Name	1 hour		3 hours		25 hours
Gene ID	Name	LFC	FDR	LFC	FDR	LFC	FDR
107848320	arogenate dehydratase/prephenate dehydratase 6, chloroplastic-like	0.63	0.51	1.69	0.00	0.35	0.79
107848097	agmatine coumaroyltransferase-2-like	1.57	0.26	2.36	0.03	1.09	0.54
107864266	chalcone synthase 1B	2.23	0.00	2.32	0.00	2.51	0.00
107871256	CHS; chalcone synthase 2	2.11	0.00	2.32	0.00	3.01	0.00
107850996	chalcone synthase J-like	2.71	0.00	3.12	0.00	2.36	0.00
107855506	dihydroflavonol-4-reductase-like	3.45	0.00	3.66	0.00	3.91	0.00
107868281	Psy; bifunctional 15-cis-phytoene synthase, chromoplastic	1.63	0.02	2.22	0.00	0.94	0.35
107867263	UPA17; growth-regulating factor 1-like	1.00	0.08	1.57	0.00	1.13	0.03
107873461	phosphomethylpyrimidine synthase, chloroplastic	1.64	0.00	2.45	0.00	1.75	0.00
107847937	pyruvate dehydrogenase E1 component subunit beta-1, mitochondrial-like	1.67	0.00	1.59	0.01	1.34	0.03
107877344	protein STRICTOSIDINE SYNTHASE-LIKE 10-like	-2.41	0.00	-1.66	0.00	-2.26	0.00
107875470	probable pyridoxal 5'-phosphate synthase subunit PDX1	1.08	0.00	1.59	0.00	0.81	0.05
107859942	adenylosuccinate synthetase 2, chloroplastic	-1.99	0.00	-1.65	0.00	-2.82	0.00
107841181	beta-amyrin synthase-like	-1.53	0.00	-1.85	0.00	-2.11	0.00
107850683	inositol-3-phosphate synthase	2.64	0.00	1.95	0.00	1.96	0.00
107864208	phosphate transporter PHO1	1.72	0.00	1.51	0.01	1.53	0.01
107863162	RNA polymerase sigma factor sigE, chloroplastic/mitochondrial	1.50	0.00	1.90	0.00	1.11	0.02
107873218	probable methionine--tRNA ligase	-1.67	0.05	-1.98	0.02	-1.10	0.27

LFC: log2 fold changes

regulated genes, the MAPK signaling pathway only was statistically significant (Fig. 6B). Moreover, the enrichment of ten categories was observed at 25 h for up-regulated genes, in which flavonoids biosynthesis and circadian rhythm-plant were significant (Fig. 7A); in contrast, 10 categories were found for down-regulated genes, but only sesquiterpenoid and triterpenoid biosynthesis were statistically significant (Fig. 7B).

Discussion

In this study, we analyzed the gene expression profile in response to combined UV-B and cold at 1, 3 and 25 h after stress exposure. The GO enrichment allowed to classify DEG into categories related to hormones, ROS oxidation, pathogens, plasma membranes and compound transport, cell wall and chloroplasts. We identified in response to stress category, three cytochrome P450 genes were up-regulated under combined stress at all times, which have been found associated to the regulation of hormone biosynthesis such as abscisic acid. These results may suggest that abscisic acid signaling leads to the maintenance of the photosynthetic activity, antioxidant enzymes activation and osmoprotectant accumulation during the combined stress of UV-B and cold (Peleg and Blumwald 2011). Moreover, we found genes associated to regulate ROS oxidation, such as catalase and peroxidase, which have been reported to be up-regulated in C. annuum subjected to cold showing protective activity (Ou et al. 2015). Interestingly, down-regulated genes (107871378, 107872419, 107875683, 107841124, 107859803, 107859806, 107856465, 107860257 and 107861184) classified within hydrolase activity were observed, they have been reported in response to pathogens, while cold triggers a negative interaction pathogen-

Fig. 4: GO enrichment analysis of genes differentially expressed at 1, 3 and 25 h in response to combined stress of UV-B radiation and cold. The categories with a (*) are statistically significant (α ≦ 0.05)

Fig. 5: Analysis of the differential genes at 1 h by KEGG enrichment map. A) Up-regulated genes, (B) down-regulated genes. The x-axis indicates the enrichment factor, and the y-axis shows the KEGG pathway. The colour of the dot represents the adjusted P - value and the size of the dot represents the number of genes

Fig. 6: Analysis of the differential genes at 3 h by KEGG enrichment map. A) Up-regulated genes B) down-regulated genes. The x-axis indicates the enrichment factor, and the y-axis shows the KEGG pathway. The colour of the dot represents the adjusted P - value and the size of the dot represents the number of genes

Fig. 7: Analysis of the differential genes at 25 h by KEGG enrichment map. A) Up-regulated genes B) down-regulated genes. The x-axis indicates the enrichment factor, and the y-axis shows the KEGG pathway. The colour of the dot represents the adjusted P - value and the size of the dot represents the number of genes

defense pathways, and UV-B radiation has been described to promote a positive interaction (Du et al. 2011; Fan et al. 2015), which may suggest that the combination of UV-B radiation and cold significantly altered the signaling networks related to pathogens, leading to the suppression of defense responses and increasing plant stem susceptibility. In addition, six genes related to plasma membranes and transport of compounds were found, where two ABC transporters were up-regulated, and involved in the transport of phytohormones, heavy metals, lipids, chlorophyll catabolites, secondary metabolites and xenobiotics (Nagy et al. 2009). Likewise, the up-expression of dynein light chain indicated activity associated with the cell membrane, acting as kinesins that transport proteins through the microtubules, from the membrane to the nucleus or vice versa (Li et al. 2018). These results indicate that compound transporter genes alleviate the disruption of osmotic and ionic homeostasis caused by UV-B and cold radiation. And up-regulation of phospholipase D and lipid phosphate phosphatase genes were observed, the phospholipase D is associated with the hydrolysis of membrane lipids and the increase of phosphatidic acid (PA) content (Li et al. 2004), and lipid phosphate phosphatase gene transforms substrates such as diacylglycerol pyrophosphate to PA and PA to diacylglycerol (Pierrugues et al. 2001). The increase in the expression of these genes at 1 and 3 h suggests high activity in the signaling of UV-B radiation and cold.

Studies have demonstrated that plants under various stresses (cold, drought, flooding and radiation) generate changes in the turgor, expansion, flexibility and rigidity of cell wall (Sasidharan et al. 2011). In this study, we detected ten down-regulated genes (107840985, 107854898, 107879143, 107840962, 107850683, 107840943, 107867324, 107878490, 107860149, 107847799), which participates in the modification and reconstruction of the cell wall, using xylan, arabinoxylan, arabinose and 1,3-β-Glucan as a substrate (Oono et al. 2006; Reboul et al. 2011). These findings indicate that the development of the stems is largely modulated by genes identified in the carbohydrate metabolic process, also it has been observed that the down-regulation of these genes limits development in pea (Lucau-Danila et al. 2012). Finally, we identified genes related to protection of chloroplasts, photoreceptor activity and flavonoid biosynthesis within response to abiotic stimulus and biosynthetic process categories. We found one sigma factor gene that was up-regulated, which regulates the transcription of chloroplast genes for the core proteins of photosystem II (Hanaoka et al. 2012); and four dehydrins, that regulate the relative loss of electrolytes, production of reactive oxygen species and chlorophyll content (Zhang et al. 2020). Moreover, 5 genes with photoreceptor activity were up-regulated, such as one stress enhanced gene that is early activated upon UV-B radiation exposure playing a photoprotective role in the thylakoid membrane (Mackerness et al. 1999), adagio-3 gene related to a photoreceptor activity to measure the duration of the day (photoperiod) (Imaizumi et al. 2003) and three UVR8 receptors, that control transcriptional responses induced by UV-B radiation (Vandenbussche et al. 2014). These findings suggest that there is an early perception of UV-B radiation at 3 h after combined stress exposure.

We examined the biochemical metabolic pathways that were affected by differential genes by KEGG enrichment analysis. Based on our results, we observed that most of up-regulated genes grouped into the flavonoids biosynthesis and circadian rhythm-plant at all sampling times. In C. annuum, an increased content of flavonoids has been observed in response to the combination by UV-B radiation and cold, maybe participating as antioxidant and UV-B absorbing compounds (León-Chan et al. 2017). We found that over time gene up-regulation was maintained in relation to products such as pinocembrin chalcone, phloretin, naringenin chalcone, 7,4’-dihydroxyflavone, apigenin and luteolin. While only in the 3 h treatment, genes related to caffeoyl-CoA were present in the production of lignin and intermediate of luteolin biosynthesis were up-regulated. At 25 h, up-expression of genes related to metabolites such as galangin, fustin, kaempferol, quercetin and myricetin were observed, which indicates that the synthesis of various flavonoids could be crucial for the protection of the plant during the first 25 h of stress. Circadian rhythm-plant was also observed at all times, Duan et al. (2014) reported that in rice abiotic stress response pathways altered the circadian clock. Interestingly, the 3 h treatment presented the over-expression of COP1 and FKF1, FKF1 works as a photoperiodic receptor for blue light (Imaizumi et al. 2003), while COP1 imports UVR8 to the nucleus from the cytosol (Yin et al. 2016), which is a UV-B specific signaling component that binds to chromatin through histones and regulates UV protection by organizing expression of a variety of genes (Rizzini et al. 2011). On the other hand, the inhibited genes FLS2, MKK9, CHIB and PYL were enriched the MAPK signaling pathway at 3 h, where FLS2 participates in the stomatal closure, a mechanism used to reduce bacterial entry into plant tissues (Mersmann et al. 2010). The MKK9 gene is related to cell death and delayed senescence in the leaves in Arabidopsis (Zhou et al. 2009). The CHIB gene has been observed in leaves and stems of sweet pepper after being infected with X. campestris pv. vesicatoria and Phytophtora capsici (Hong et al. 2000). This suggests that at 3 h after treatment the pepper plants show greater sensitivity to infection by pathogens.

Conclusion

We performed the transcriptomic analysis of the combined effect of UV-B radiation and cold on stems of C. annuum L. after stress exposure at 1, 3 and 25 h. We identified the induction of genes related to abscisic acid biosynthesis at 1 h. Furthermore, we can infer that after 3 h there is the greatest susceptibility to pathogens. We also observed that in response to combined stress, genes associated to flavonoid biosynthesis are induced at 1 h after treatment. These data will be very useful genetic resource to analyze the resistance of peppers to cold and UV-B radiation. Furthermore, further studies are needed to confirm the roles of the candidate genes in the identified processes.

Acknowledgements

This work was supported by FOSEC SEP-INVESTIGACIÓN BÁSICA, Proyecto No. A1-S-8466. Cátedras CONACYT: Proyecto No. 784. Lightbourn Research. Convenio: 589683, Proyecto: Análisis Transcripcional de Pimiento Morrón (Capsicum annuum L.) bajo estrés abiótico. The authors thank Q.F.B. Jesús Héctor Carrillo Yáñez for critical technical assistance.

Author Contributions

JLF, BH and LLR conceived, designed and coordinated the study. BMM and RLC carried out the experimentation. CV, JLF, APT, OCR and HMM analyzed the results. Contributed reagents/materials/analysis tools: JLF, LLR and BH. CV and BH edited the English grammar of the manuscript. All authors wrote, read and approved the final manuscript.

Conflicts of Interest

All other authors declare no conflicts of interest.

Data Availability

Data presented in this study are available on fair request to the corresponding author.

Ethics Approval

Not applicable.

References

Airaki M, M Leterrier, R Mateos, R Valderrama, M Chaki, J Barroso, LD Río, J Palma, F Corpas (2011). Metabolism of reactive oxygen species and reactive nitrogen species in pepper (Capsicum annuum L.) plants under low temperature stress. Plant Cell Environ 35:281‒295

Anders S, W Huber (2010). Differential Expression for Sequence Count Data. Genom Biol 11:4310‒1415

Atkinson N, P Urwin (2012). The interaction of plant biotic and abiotic stresses: From genes to the field. J Exp Bot 63:3523‒3543

Chugh S, S Sharma, A Rustagi, P Kumari, A Agrawal, D Kumar (2018). Enhancing cold tolerance in horticultural plants using in vitro approaches. In: Abiotic Stress-Mediated Sensing and Signaling in Plants: An Omics Perspective, 1^st Edition, pp:225-241. Springer, Singapore

Du H, Y Liang, K Pei, K Ma (2011). UV radiation-responsive proteins in rice leaves: A proteomic analysis. Plant Cell Physiol 52:306‒316

Duan M, P Huang, X Yuan, H Chen, J Huang, H Zhang (2014). CMYB1 Encoding a MYB transcriptional activator is involved in abiotic stress and circadian rhythm in rice. Sci World J 2014; Article 178038

Fan S, L Jiang, J Wu, L Dong, Q Cheng, P Xu, S Zhang (2015). A novel pathogenesis-related class 10 protein Gly m 4l, increases resistance upon Phytophthora sojae infection in soybean (Glycine max [L.] Merr.). PLoS One 10; Article e0140364

Guo WL, RG Chen, ZH Gong, YX Yin, SS Ahmed, YM He (2012). exogenous abscisic acid increases antioxidant enzymes and related gene expression in pepper (Capsicum Annuum) leaves subjected to chilling stress. Gene Mol Res 11:63‒80

Hanaoka M, M Kato, M Anma, K Tanaka (2012). SIG1, a sigma factor for the chloroplast RNA polymerase, differently associates with multiple DNA regions in the chloroplast chromosomes in vivo. Intl J Mol Sci 13:12182‒12194

Hong J, H Jung, Y Kim, B Hwang (2000). Pepper gene encoding a basic class II chitinase is inducible by pathogen and ethephon. Plant Sci J 159:39‒49

Hou X, H Zhang, S Liu, X Wang, Y Zhang, Y Meng, D Luo, R Chen (2020). The NAC transcription factor CaNAC064 is a regulator of cold stress tolerance in peppers. Plant Sci 291; Article 110346

Hwang E, K Kim, S Park, M Jeong, M Byun, H Kwon (2005). Expression profiles of hot pepper (capsicum annuum) genes under cold stress conditions. J Biosci 30:657‒667

Imaizumi T, H Tran, T Swartz, W Briggs, S Kay (2003). FKF1 is essential for photoperiodic-specific light signalling in Arabidopsis. Nature 426:302‒306

Koç E, C Işlek, AS Üstün (2010). Effect of cold on protein, proline, phenolic compounds and chlorophyll content of two pepper (Capsicum annuum L.) varieties. Gazi Univ J Sci 23:1‒6

Lai Y, B Xu, L He, M Lin, L Cao, S Mou, S He (2011). Differential gene expression in pepper (Capsicum annuum) exposed to UV-B. Ind J Exp Biol 49:429‒437

León-Chan RG, LA Lightbourn-Rojas, M López-Meyer, L Amarillas, J Basilio, TF Martínez-bastidas, C Villicaña, J León-Félix (2020). Differential gene expression of anthocyanin biosynthetic genes under low temperature and ultraviolet-B radiation in bell pepper (Capsicum annuum). Intl J Agric Biol 23:501‒508

León-Chan RG, M López-Meyer, T Osuna-Enciso, J Sañudo-Barajas, J Heredia, J León-Félix (2017). Low temperature and ultraviolet-B radiation affect chlorophyll content and induce the accumulation of UV-B-absorbing and antioxidant compounds in bell pepper (Capsicum anuum) plants. Environ Exp Bot 139:143‒151

Li J, D Yu, G Qanmber, L Lu, L Wang, L Zheng, Z Liu, H Wu, X Liu, Q Chen, F Li, Z Yang (2018). GhKLCR1, a kinesin light chain-related gene, induces drought-stress sensitivity in Arabidopsis. Sci Chin Life Sci 62:63‒75

Li W, M Li, W Zhang, R Welti, X Wang (2004). The plasma membrane-bound phospholipase Dδ enhances freezing tolerance in Arabidopsis thaliana. Nat Biotechnol 22:427‒433

Love MI, W Huber, S Anders (2014). Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol 15:1‒21

Lucau-Danila A, C Toitot, E Goulas, AS Blervacq, D Hot, N Bahrman, H Sellier, I Lejeune-Hénaut, B Delbreil (2012). Transcriptome analysis in pea allows to distinguish chilling and acclimation mechanisms. Plant Physiol Biochnol 58:236–244

Mackerness SA, SL Surplus, P Blake, CF John, V Buchanan-Wollaston, BR Jordan, B Thomas (1999). ultraviolet-b-induced stress and changes in gene expression in arabidopsis thaliana: role of signalling pathways controlled by jasmonic acid, ethylene and reactive oxygen species. Plant Cell Environ 22:1413‒1423

Mahdavian K, M Ghorbanli, K Kalantari (2008). The effects of ultraviolet radiation on the contents of chlorophyll, flavonoid, anthocyanin and proline in Capsicum annuum L. Turk J Bot 32:25‒33

Mercado J, M Reid, V Valpuesta, M Quesada (1997). Metabolic changes and susceptibility to chilling stress in capsicum annuum plants grown at suboptimal temperature. Funct Plant Biol 24:759‒767

Mersmann S, G Bourdais, S Rietz, S Robatzek (2010). ethylene signaling regulates accumulation of the fls2 receptor and is required for the oxidative burst contributing to plant immunity. Plant Physiol 154:391‒400

Mittler R, E Blumwald (2010). Genetic engineering for modern agriculture: challenges and perspectives. Annu Rev Plant Biol 61:443‒462

Nagy R, H Grob, B Weder, P Green, M Klein, A Frelet-Barrand, J Schjoerring, C Brearley, E Martinoia (2009). The Arabidopsis ATP-binding cassette protein AtMRP5/AtABCC5 is a high affinity inositol hexakisphosphate transporter involved in guard cell signaling and phytate storage. J Biol Chem 284:33614‒33622

Oono Y, M Seki, M Satou, K Iida, K Akiyama, T Sakurai, M Fujita, K Yamaguchi-Shinozaki, K Shinozaki (2006). Monitoring expression profiles of Arabidopsis genes during cold acclimation and deacclimation using DNA microarrays. Funct Integr Genomics 6:212‒234

Ou L, G Wei, Z Zhang, X Dai, X Zou (2015). Effects of low temperature and low irradiance on the physiological characteristics and related gene expression of different pepper species. Photosynthetica 53:85‒94

Peleg Z, E Blumwald (2011). Hormone balance and abiotic stress tolerance in crop plants. Curr Opin Plant Biol 14:290‒295

Pierrugues O, C Brutesco, J Oshiro, M Gouy, Y Deveaux, GM Carman, P Thuriaux, M Kazmaier, (2001). Lipid phosphate phosphatases in Arabidopsis regulation of the AtLPP1 rene in response to stress. J Biol Chem 276:20300‒20308

Pressman E, R Shaked, N Firon (2006). Exposing pepper plants to high day temperatures prevents the adverse low night temperature symptoms. Physiol Plantarum 126:618‒626

Reboul R, C Geserick, M Pabst, B Frey, D Wittmann, U Lütz-Meindl, R Léonard, R Tenhaken (2011). Down-regulation of UDP-glucuronic acid biosynthesis leads to swollen plant cell walls and severe developmental defects associated with changes in pectic polysaccharides. J Biol Chem 286:39982‒39992

Rizzini L, J Favory, C Cloix, D Faggionato, A O'Hara, E Kaiserli, R Baumeister, E Schafer, F Nagy, G Jenkins, R Ulm (2011). Perception of UV-B by the Arabidopsis UVR8 protein. Science 332:103‒106

Rodríguez-Calzada T, M Qian, Å Strid, S Neugart, M Schreiner, I Torres-Pacheco, R Guevara-González (2019). Effect of UV-B radiation on morphology, phenolic compound production, gene expression, and subsequent drought stress responses in chili pepper (Capsicum annuum L.). Plant Physiol Biochem 134:94‒102

Sasidharan R, L Voesenek, R Pierik (2011). Cell wall modifying proteins mediate plant acclimatization to biotic and abiotic stresses. Crit Rev Plant Sci 30:548‒562

Suzuki N, R Rivero, V Shulaev, E Blumwald, R Mittler (2014). Abiotic and biotic stress combinations. New Phytol 203:32‒43

Vandenbussche F, K Tilbrook, AC Fierro, K Marchal, D Poelman, DVD Straeten, R Ulm (2014). Photoreceptor-mediated bending towards UV-B in Arabidopsis. Mol Plant 7:1041‒1052

Yi S, J Kim, Y Joung, S Lee, W Kim, S Yu, D Choi (2004). The pepper transcription factor capf1 confers pathogen and freezing tolerance in Arabidopsis. Plant Physiol 136:2862‒2874

Yin R, M Skvortsova, S Loubéry, R Ulm (2016). COP1 is required for UV-B–induced nuclear accumulation of the UVR8 photoreceptor. Proc Natl Acad Sci USA 113:4415‒4422

Zhang HF, SY Liu, JH Ma, XK Wang, SU Haq, YC Meng, YM Zhang, RG Chen (2020). CaDHN4, a salt and cold stress-responsive dehydrin gene from pepper decreases abscisic acid sensitivity in Arabidopsis. Intl J Mol Sci 2; Article 26

Zhou C, Z Cai, Y Guo, S Gan (2009). An Arabidopsis mitogen-activated protein kinase cascade, mkk9-mpk6, plays a role in leaf senescence. Plant Physiol 150:167‒177